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Abstract 

This paper describes the results of our theoretical and numerical studies of hy- 
drodynamic interactions in a suspension of spherical particles confined between two 
parallel planar walls, under creeping-flow conditions. We propose a novel algorithm 
for accurate evaluation of the many-particle friction matrix in this system — no such 
algorithm has been available so far. 

Our approach involves expanding the fluid velocity field into spherical and Carte- 
sian fundamental sets of Stokes flows. The interaction of the fluid with the particles 
is described using the spherical basis fields; the flow scattered with the walls is ex- 
pressed in terms of the Cartesian fundamental solutions. At the core of our method 
are transformation relations between the spherical and Cartesian basis sets. These 
transformations allow us to describe the flow field in a system that involves both 
the walls and particles. 

We used our accurate numerical results to test the single-wall superposition ap- 
proximation for the hydrodynamic friction matrix. The approximation yields fair 
results for quantities dominated by single particle contributions, but it fails to de- 
scribe collective phenomena, such as a large transverse resistance coefficient for 
linear arrays of spheres. 



1 Introduction 



Equilibrium and nonequilibrium behavior of colloidal suspensions in confined 
geometries has recently been extensively discussed. Examples of recent papers 
include experimental studies of particle deposition on chemically patterned 
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planar walls 1], investigations of collective dynamics in quasi-bidimensional 
suspensions in slit pores 0, 0, 0], and observations of drainage behavior of 
particle- stabilized thin liquid films 0. The research has been stimulated, in 
part, by emerging applications — such as microfluidic devices and production 
of photonic materials by self-assembly of colloidal crystals jgl Q . The investi- 
gations have also been considerably influenced by development of new exper- 
imental techniques, including the evanescent-wave microscopy S jl, comput- 
erized video microscopy 1^, HI, 0, E3 5 an d optical tweezers |l3| . 



While the equilibrium structure of confined colloidal suspensions is fully de- 
termined by the particle-wall and interparticle interaction potentials, the dy- 
namics is also significantly affected by the many-body hydrodynamic forces. 
The effect of the hydrodynamic interactions on particle motion can be ex- 



pressed in terms of the iV-particle friction and mobility matrices [1^, which 
depend on the particle positions and the wall geometry. 

For spherical particles in an unbounded space, efficient algorithms for evalu- 
ation of the friction and mobility matrices have been developed |l5[ [lH 17 . 

The algorithms combine multipolar expansion methods with the lu- 
brication approximation for particles in close proximity [15(. This approach 
has recently been generalized by Cichocki et al. 2(J 21 1 to systems of particles 



bounded by a single planar wall; the particle-wall hydrodynamic interactions 
were included using the image representation of the flow reflected from the 
wall H- 



Much less progress has been made for suspensions confined between two planar 
walls (e.g, in a slit pore, or between two glass plates). For a single particle, 
several ad hoc approximations for the mobility matrix have been proposed 
[22I l23l |. and numerical results obtained by boundary-integral methods are 
available (24 , 25[ 26| . Recently, we have developed an exact image representa- 



tion of the flow between two walls [27j , which allows accurate evaluation of the 
single-particle friction matrix by a multipolar expansion technique. However, 
none of the above methods has been generalized to multiparticle systems, 
due to a large numerical cost of boundary-integral calculations or the slow 
convergence of the image solutions. 

Two extensions of the free-space Stokesian-dynamics algorithm to wall 
bounded systems have been proposed by Brady and his collaborators (H, 
3oj | . In the first approach the walls are discretized [2^], and in the second 



they are modeled as static, closely packed arrays of spheres 2^, 3(j . The first 
method has not been further explored. The results obtained using the second 
method are only qualitative, because the walls are porous and rough. 

To overcome the above-mentioned problems, we adopt here an alternative ap- 
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proach, based on Fourier analysis of the flow field reflected from the walls. 2 
According to our method, the flow field in the system is expanded using two 
basis sets of solutions of Stokes equations — the spherical and the Cartesian 
basis. The spherical basis is used for a description of the flow field scattered 
from the particles, and the Cartesian basis is appropriate for the wall geometry. 
The key result of our study is a set of transformation formulas for conversion 
between the spherical and Cartesian representations. The transformation for- 
mulas allow to evaluate the spherical matrix elements of the Green function 
for Stokes flow in the presence of the walls in terms of simple two-dimensional 
Fourier integrals. 

The results of our theoretical analysis have been implemented in a numer- 
ical procedure for evaluating multi-particle hydrodynamic interactions in a 
suspension of spheres confined between two planar walls. The procedure com- 
bines the expansions of the flow field into the spherical and Cartesian basis 
fields with the two-particle superposition approximation for the friction ma- 
trix, in order to include slowly convergent lubrication corrections. Since the 
force multipoles induced on particle surfaces are included to arbitrary order, 
highly accurate results are obtained. 

Examples of numerical results for two-particle and many-particle systems are 
provided. In particular, our results illustrate the role of the far-field flow pro- 
duced in the space between the walls by the moving particles. We show that 
the single-wall superposition approximation does not correctly describe the 
far-field flow, and thus it fails to capture some important collective phenomena 
such as the increased hydrodynamic resistance due to the backflow produced 
by the moving particles. 

This paper is organized as follows. In Section 2 we summarize the induced- 
force formulation of the Stokes-flow equations for a multiparticle system in 
the wall presence. In Section 3 the induced-force equations are transformed 
into an infinite array of algebraic equations for the induced-force multipoles, 
using a multipolar expansion of Stokes flow. 

Our main theoretical results are presented in Sections 4-7. In Section 4 the 
Cartesian basis sets of Stokes flows are defined, and the transformation formu- 
las for conversion between the Cartesian and the spherical multipolar bases 
sets are derived. The displacement and conversion formulas are then used 
to obtain two-dimensional Fourier representations of the matrix elements of 
Green operator for infinite space (Section 5), halfspace bounded by a single 
wall (Section 6), and a region bounded by two parallel planar walls (Section 
7). 



2 Recently, investigations along similar lines have also been reported by Jones 31 
0. 
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A numerical algorithm for computation of hydro dynamic interactions in a 
suspension of spheres confined in a region bounded by two parallel walls is 
described in Section 8. Numerical examples of the friction matrix, evaluated 
using this algorithm, are given in Section 9. Directions for further development 
of our method are indicated in the concluding Section 10. Some technical 
details are presented in Appendices A-D. 



2 Induced-force formulation 



We consider a suspension of N spherical particles of radius a moving in an 
incompressible Newtonian fluid of viscosity rj. The suspension is bounded by 
a single planar wall or two parallel planar walls. The creeping-flow condi- 
tions are assumed; therefore, the fluid flow in the system depends only on 
the instantaneous particle configuration and velocities. The configuration is 
described by the positions (Ri, . . . , Rjv) of particle centers. The translational 
and rotational velocities of the particles are Uj and Hi, where i — 1, . . . , N. 

The effect of the suspended particles on the surrounding fluid can be described 
in terms of the induced force distributions on the particle surfaces 

F l (r)=a- 2 5(r l -a)f l (r), (1) 

where 

r, = r - R (2) 
and r, = |rj|. By definition of the induced force, the flow field 

N 

v(r) = v cxt T(r, r') • F 4 (r') dr' (3) 

i=i 

is identical to the velocity field in the presence of the particles 33, 34, 35l . 
In the above equation, v ext denotes the imposed flow, and the integral term 
describes the flow generated by the induced forces. Here 

T(r,r / ) = To(r-r , ) + T / (r,r') (4) 

is the Green function for the Stokes flow in the presence of the boundaries, 

T (r) = -^(I + ff) (5) 

denotes the Oseen tensor (where i is the identity tensor, and r = r/r), and 
T'(r, r') describes the flow reflected from the walls. 

The induced force distribution Fj on the surface of particle i and the flow v™ 
incident to this particle are linearly related. The relation can be expressed in 
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the form 

F^-Z^vf-vf), (6) 

where 

vf(r) = V l + n i xv l (7) 
denotes the rigid-body velocity field corresponding to the particle motion, and 

vj n = v? - vf (8) 
is the incident flow in the reference frame moving with the particle. The Stokes 
flow field (8) is fully determined by its boundary value on the particle surface 
Si and the condition that vf 1 is nonsingular in the region occupied by the 
particle. Thus, (6) can be interpreted as a linear functional relation between 
the vector fields (1) and (8) specified on the surface Si. Since a nonzero flow (8) 
always produces a nonzero force distribution Fj, relation (6) can be inverted 

v in_ v rb = _ [Zr l Fi](r)) reS ._ (Q) 

For specific particle models, explicit expressions for the operator Zj are ob- 
tained by solving Stokes equations for an isolated particle subject to an ex- 
ternal flow in an unbounded fluid 36l 37. iil . 

The flow yf 1 incident to a particle i in a multiparticle system is defined by the 
equation 

v(r) =<(!•)+<>), (10) 
where v(r) is the total flow (3), and 



v out 



» = J Totr-rO-F^rOdr' (11) 



represents the flow scattered by the considered particle. By collecting relations 
(9)— (11) we obtain the expression 

v(r)=vf(r)-[Z7 1 F,](r) + J T (r - r') • F,(r') dr', r G Si, (12) 

for the flow at the surface Si of the particle i. For rigid spheres, the velocity 
field v(r) in equation (12) satisfies the no-slip boundary condition 

v(r)=vf(r), reS,. (13) 

Accordingly, we have the identity 21 



[Zr 1 F l ](r) = J T (r-r')-F ?; (r')dr', r G Si (14) 
for such particles. 

By combining expressions (3) and (12), we get the boundary-integral equation 
for the induced force densities Fj, 
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[Zr l P*](r) +E /K 1 - %) T o(r - r) + T'(r, r')] • F,(r') dr' = vf> (r) - v<* 

r G 

In the following sections, equation (15) is transformed into an infinite set 
of algebraic equations for the multipole moments of the induced force, with 
coefficients expressed in terms of two-dimensional Fourier integrals. 



3 Matrix representation 



3.1 Spherical basis 



The matrix representation of equation (15) is obtained by expanding fluid 
velocity fields into sets of fundamental solutions of Stokes equations in spheri- 
cal coordinates, and expressing the induced-force distributions in terms of the 
corresponding force multipoles. In our analysis we employ sets of basis fields 
that are closely related to the sets introduced by Cichocki et al. |37| ; we use, 
however, a different normalization to emphasize important symmetries of the 
problem. 

The singular and nonsingular basis sets of solutions of Stokes equations v^ CT (r) 
and v z { ncr (r) (where I = 1, 2, . . .; m = —I, and er = 0, 1, 2) are defined by 
the following conditions: (i) the basis velocity fields are homogeneous functions 
of the radial variable r, 

vr mff (r) = V imff (^)r^, (16) 



vL,(r)=V^(My +CT -\ (17) 
where (r, 6, <fi) represent the vector r in spherical coordinates; (ii) the coeffi- 
cients V z ^ no .(^, (f>) and V^^, 0) are combinations of vector spherical harmon- 
ics with angular order / and azimuthal order m; and (in) the basis velocity 
fields v^ CT (r) satisfy the following hierarchies of curl relations 



and 



V /ml = ~ iV X V ImO. 

V /m2 = ~ iV X V *ml> 
V Zml = IV X V^ 2 , 



(18a) 

(18b) 
(19a) 
(19b) 
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The above identities imply that only the solutions v^ and v^ 2 have nonzero 
corresponding pressure fields, and that the solutions v^ 2 and v^ represent 
potential flows, i.e., 



Vxv 



lm2 



0. 



Vxv 



ImO 



0. 



(20a, b) 



The proportionality coefficient in the curl relations (18) and (19) is determined 
by the requirement that the ex pan sion of the Oseen tensor in basis functions 
(16) and (17) has the form 37, 39l 



^To(r-r') 



Ima 



Ima 



rv 



Vv+ ( 

^ Ima 



Ima \ 



Ima \ 



ty ~~~~~ ;r ry 



ry < ^~^ ry 



(21a, b) 



The conditions (16)-(21) determine the basis fields v^ up to a single nor- 
malization constant, which is set by an additional requirement that 



(221 



where Yi m is the normalized scalar spherical harmonic (as defined by Edmonds 
0). 

The flow fields (16) and (17) form complete sets of singular and non-singular 
solutions of Stokes equations in the representation appropriate for spherical 
symmetry. However, they do not form orthonormal sets with respect to the 
natural functional scalar product for vector fields A and B on the spherical 
surface r = b. Following the approach of Cichocki et al. 37J we thus intro- 
duce the reciprocal basis fields w^, which are defined by the orthogonality 
relations 

( S b W La I V /w) = Mmm'<W (23) 

for all values of parameter b > 0, where 



^(r) = 6-^(r-6), 



(24) 



and 



(A | B) = J A*(r) • B(r) dr. (25) 

The functions w^ CT have a similar structure to the functions v^ CT , i.e., they 
have a separable form 



w 



Ima \ 



w 



lma\ u i ' 



(26) 



W lma{ r ) 



-(l+a-1) 



(27) 
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with the coefficients W lma (9, 0) given by combinations of vector spherical 
harmonics with angular order / and azimuthal order m. Explicit relations for 
the functions V^ CT (0, 0) and W^(0, <p) in equations (16), (17), (26), and (27) 
are listed in Appendix A. 

3.2 Equations for induced-force multipole moments 

In the multipolar-representation method, the boundary- integral equation (15) 
is transformed into an infinite set of linear algebraic equations for the multipo- 
lar moments of the induced-force distributions (1). The multipolar expansion 
of the distribution F; is defined by the relation 

F 4 (r) = £ /^m^ a s (r>+ >,)■ (28) 

Ima 

The corresponding multipolar moments are given by 

/iCZma) =|viJ^(r i ).F i (r)dr 

= (<M |P*>, (29) 

consistent with the orthogonality relation (23). In the above equation we in- 
troduce the standard bra-ket notation, with an additional convention that 
|A) represents the vector field A(r) and |A(i)) denotes A(i-j). 

The linear algebraic equations for the multipolar moments of the induced 
force (29) are obtained by projecting the linear operators in the boundary- 
integral equation (15) onto the reciprocal basis (27). The resulting matrix 
representation of equation (15) can be written in the form 

N 

]T M^lma | I'm' a') fj (I'm' a') = a(lma), (30) 

j=l I'm 1 a' 

where 

Minima | I'm'a') = Z'^lma \ I'm 1 a') + G% (Ima | I'm'a') +G' ij (lma | I'm'a'), 

(31) 

and 

Z^(lma | I'm'a') = M^Wwj^i) I Zj 1 | 6 s a (j)w+ m , a ,(j)), (32) 
G%(lma | I'm'a') = (1 - <%) <5 a s «w+ mCT « | T | 6 s a (j)w+ m , a/ (j)), (33) 
G' l3 (lma | I'm'a') = (<5 a s (i)w+ a {i) | T' | if 0>? m v0')>- (34) 
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In equations (32)-(34) 



[Tb] (r) = J T(r, r') • b(r') dr', (35) 

and the bra-ket notation introduced in equation (29) is used. The matrix el- 
ements (33) and (34) are independent of the particle radius a, because the 
orthogonality relation (23) holds for all values of the parameter b. The coeffi- 
cients Ciilma) on the right side of equation (30) are defined by the expansion 

vf (r) - v cxt (r) = £ Q(/ma)v+ (36) 

Imcr 

of the imposed flow field relative to the rigid-body particle motion (7) into the 
basis functions (17) centered at the position of particle %. Inserting the above 
expression into the orthogonality relation (23) yields 

c t {lma) = (5 s a (i)w+ ma (i) | vf - v cxt ). (37) 



For a system of identical particles, the matrix elements (32) of the one-particle 
operator Z" 1 are independent of the particle label i. Since the particles are as- 
sumed to be spherical, the matrix elements (32) are diagonal in the multipolar 
orders I and m, and independent of m, 

Zr^(lmo | I'm'o') = 5 ij 5 w 5 mm iZ^ 1 (l] a \ a'). (38) 

By specifying equation (30) for a single isolated particle i in an unbounded 
fluid and using the diagonality relation (38) we obtain the linear formula 

£ Zr\l- a | a') Mima') = Ci (lma). (39) 

cr' 

Inserting the Oseen tensor in the form (21a) into equation (11) and using the 
definition (29) we also get 

v° ut (r) = tT 1 E Mrrur^iTi). (40) 

Imcr 



According to the above relation, the multipolar moments fi(lma) can be inter- 
preted as the expansion coefficient of the flow field v° ut scattered by the parti- 
cle i into the basis velocity fields (16). It follows that the matrix ^ i ~ 1 (/; o \ cr') 
relates the expansion coefficients of the incident and the scattered flows. For 
hard spheres, porous particles, spherical viscous drops, and spherical drops 
covered by an incompressible surfactant layer, explicit expressions for the ma- 



trix elements (38) are known |36j, 1371 138|. (Note, however, a different normal- 



ization of the basis functions here and in the above references, as discussed in 
Appendix A) 
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The matrix elements of the free-space Oseen operator (33) are also known, 
since they are simply linked to the elements of the displacement matrix S~* 



that was evaluated by Felderhof and Jones [4l( . To show this relation we insert 
expression (21a), specified for T (rj — r^) with r j = r — Hj and r'- = r' — PL,-, 
into equation (33), and use the orthogonality condition (23) for the fields 
centered at the position of the particle j. As the result we find 



21 



G%(lma | I'm'a') = v' 1 (S s a (z)w+ m „(*) | v7 mV ,(j)>- (41) 

The matrix element on the right side of the above equation corresponds to 
the expansion of the singular flow field v7 m ' CT / centered at the position of the 
particle j into the nonsingular basis flow fields v^ CT centered at the position 
of the particle i, 

v,w(r,) = E vL^K^Ww+^W | v t 7 mV 0")>- ( 42 ) 

Ima 

According to the definition of the displacement matrix jilj we thus have 

G%{lma | I'm'a') = rT^^R* - R;; Ima | I'm'a'). (43) 

We note that the displacement matrix S$~ , introduced above, is normalized 
differently than the matrix S + ~ defined by Felderhof and Jones 0] (cf., the 
transformation (A. 12) between the corresponding basis fields.) 

As a result of the Lorentz symmetry 

T Q (r,r') = Tt(r',r) (44) 

of the Green functions T a = T , T' (where the dagger denotes the transpose of 
the tensor) and the symmetry of the scalar product (25), the matrix elements 
(33) and (34) satisfy the reciprocal relations 



G^lma | I'm'a') = G]*{l'm'a ! \ Ima), (45) 
G'tjilma | I'm'a') = G'/^l'm'a' | Ima), (46) 

where the asterisk denotes the complex conjugate. The matrix elements (38) 
of the one-particle scattering operator have a similar symmetry, 

Zr\l-o\a') = Z.\l-a'\a). (47) 

The matrix elements (47) are real, due to the diagonality (38) of the matrix 
(32) in the azimuthal number m. 
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3.3 Matrix notation 



In what follows we will use a compact matrix notation in the three-dimensional 
linear space with the components corresponding to the indices a = 0,1,2 
that identify the tensorial character of the basis flow fields (16) and (17). 
Accordingly, the matrices with the elements (32)-(34) will be denoted by 
Z^ 1 (/m | I'm'), G?(Zm | I'm'), and G'ij(lm \ I'm'), respectively; the matrices 
with the elements Z~ l (l;a \ a') and Sg~~(Rj — Hj-Jma \ I'm'a') will be de- 
noted by Z~ X (Z) and Sg~(Rj — R,-; Im \ I'm'). A similar convention will be used 
for three-dimensional column vectors representing quantities with a single in- 
dex a (such as the induced-force multipolar amplitudes). With this notation, 
equation (30) can be written in the form 

N 

E M n ( lm I l ' m ') • f j ( l ' m> ) = °i (H > ( 48 ) 

j=l I'm' 

with the matrix M given by the relation 

Mijilm | I'm') = ^n-L-Zr'p) + G%{lm \ I'm ) + G'^lm | I'm'), (49) 

according to expressions (31) and (38). In the above equations, fj(/m) and 
Ci(lm) are column vectors with components fi(lma) and Cj(Zma), and the dot 
represents the matrix multiplication. 

An analogous matrix notation will be used in the Cartesian representation, 
which is introduced in the following section. 



4 Cartesian basis 



Using the force-multipole equations (48) to determine the hydrodynamic fric- 
tion matrix in a suspension bounded by planar walls involves evaluation of the 
spherical matrix elements (34) of the Green function T'(r, r') that describes 
the flow field in the bounded domain. For a single wall the matrix elements 
were calculated by Cichocki et al. 2(1 21] using a multipolar-image represen- 
tation of the flow reflected from a planar boundary. For a suspension confined 
between two parallel walls, the matrix elements (34) can be evaluated using 
the image representation derived by Bhattacharya and Blawzdziewicz |27f : 
however such calculations are inefficient due to convergence problems. Here 
we propose an alternative approach, in which the matrix elements (48) are 
determined by means of Cartesian representation of the flow fields, consistent 
with the wall geometry. 
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In what follows we assume that the walls are normal to the axis z in the Carte- 
sian coordinate system (x,y,z). The corresponding Cartesian unit vectors are 
denoted e x ,e y ,e z . 



4-1 Definition of Cartesian basis flow fields 

We introduce two basis sets of Stokes flows v kcr and v k(T , defined by the ex- 
pressions 

v k0 (r) = (32tt 2 )- 1/2 [i(l-2kz)k+(l + 2kz)e z ]k' 1/2 ^- p - kz , (50a) 

v kl (r) = (87r 2 )- 1 / 2 (k x e z )k~ 1 ' 2 ^' p ~ k \ (50b) 

v k2 (r) = (32tt 2 )- 1/2 (ik - e 2 )AT 1/2 e ik " p - kz , (50c) 

and 

v+(r) = (327T 2 )- 1 / 2 (ik + e z )k- 1 ' 2 e ^ p+kz , (51a) 

v+(r) = (87r 2 )- 1 / 2 (k x e^k' 1 ' 2 e ik ' ' p+kz : , (51b) 

v+ (r) = (327r 2 )- 1 / 2 [i(l + 2kz)k-(l-2kz)e 2 ] k~ 1 ' 2 ^ p+kz . (51c) 
The pressure fields corresponding to the flows (50) and (51) are 
p k0 (r) = (2n 2 )- 1 / 2 r ] k 1 / 2 e^ p - kz , p+ (r) = {2ir 2 y 1/2 ri k l/2 e ik ' p+kz , (52) 



and 

Pki( r ) = Pk 2 ( r ) = Pko( r ) = Pki( r ) = 0. 
In the above relations 

p = xe x + ye y 

is the projection of the vector r onto the x-y plane, and 

k k x &x ~t~ kyGy 



(53) 
(54) 

(55) 



is the corresponding two-dimensional wave vector. Furthermore, k = k/k and 
fc = 

The basis sets (50) and (51) are determined by the following conditions: Firstly, 
each basis flow field v ko . corresponds to a single lateral Fourier mode; secondly, 
the velocity fields v kCT (r) vanish for z — > 00 and v kcr (r) for z — > —00; thirdly, 
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the basis fields v k(7 are obtained from v k2 _ CT by reflection with respect to the 
plane x-y. The fourth condition is the set of curl relations 

v kl = -\ik~ l V x v k0 , (56a) 



'k2 



-\ik X V x v kl , 



and 



v+ = hk- 1 ^ X v+, 



(56b) 
(57a) 



v+ = \ik- l V X v+, (57b) 

by analogy to the expressions (18)-(20) for the spherical basis. Relations (56) 
and (57) imply 

Vxv k2 = 0, Vxv+=0. (58a, b) 

The final condition is the requirement that the basis fields (50) and (51) satisfy 
the identity 



' /dk£v kCT (r)v+;(r'), z>z>, 



^T (r-r') 



(59a, b) 



/dkEv k + ff (r)v k ;(r'), 

Jk 



Z < Z, 



where the integration is over the two-dimensional Fourier space (55). The 
identity (59) is analogous to the representation (21) of the Oseen tensor in the 
spherical basis. It can be verified by showing that 



vL(r)v^(r') = r ? t (k,z-/) e ik '^'), 



(60) 



where 



r)T (k,z) = ^[2i-(l + A;|2|)kk-i^(ke 2 + e 2 k)-(l-A;|2|)e 2 e 2 ]A;- 1 e-^l, 

(61) 

is the two-dimensional Fourier transform in the x-y plane of the Oseen tensor, 
to(k,^) = ^p/T (r)e- ik ^dp. (62) 

The reciprocal sets of vector fields w k(T that correspond to the Cartesian basis 
sets v k(7 are defined by the orthogonality relations 



(^w± |v^> = < 5(k-k0^, 
which hold for all values of the parameter h, where 

5 c h (v) = 5(z-h). 



(63) 



(64) 
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Inserting expressions (50) and (51) into the above relations yields 

w±(r) = 4*v£,(r), w±(r) = 2fcv£(r), w± (r) = 4*v£(r). (65) 

^.2 Displacement theorem for Cartesian basis 



The Cartesian basis fields (50) and (51) centered at different points Ri and R2 
are related by simple displacement transformations. Due to the translational 
invariance, the transformations are diagonal in the wave vector k, 



v ka( r 2)=^v ka ,(n)5 c (Ri2,k ;( j' I a), 

a' 

v^(r 2 )=E v w( r i)^c + (Ri2,k ;( 7' I a), 



(66a) 
(66b) 



where ri = r — R l5 r 2 = r — R 2 , and R i2 = Ri — R 2 . Using the orthogonality 
condition (63) and the completeness of the Cartesian basis sets we get 



<5 c (l)w± CT ,(l) I v± (2)> = 5(k-k')^ ± (R 12 ,k; ( T' I a). 



(67) 



Relations (66) and the expression (67) for the elements 5'c ± (Ri 2 , k; a' \ a) of 
the Cartesian displacement matrix Sq ± (Ri 2 , k) are analogous to the displace- 
ment formulas (41)-(43) for the spherical basis fields. 

An analysis of equations (50) and (51) indicates that the matrices Sc ± (Ri2, k) 
can be written in the factorized form 



S ±± 



where 
and 



c iRi2,k) 

Rl2 — P12 + Zl2&z: 



(68) 
(69) 



Sc~(kZ) 



1 
1 
-2kZ 1 



S^(kZ) 



1 2kZ 
1 
1 



(70) 



It is also easy to verify that the matrices obey the group property 

S±±(R + R, k) = S±±(R, k) • S^(R, k) (71) 

with 



S^(0,k) = I, 



(72) 
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and that they satisfy the symmetry relation 

S+ + (R,k) = [S c -(-R,k)]t, (73) 
where the dagger denotes the Hermitian conjugate. 

4-3 Transformations between Cartesian and spherical basis sets 



One of the key results of our study is a set of transformation relations between 
the spherical basis fields (16) and (17) and the Cartesian basis fields (50) and 
(51). We focus on the transformations 

vf m ,(r 2 ) = / dk'Ev^,( ri )(5 c (l)w^,(l) I v r _(2)) (74) 

a' 

and 

vi^(r 2 ) = £ v+ m , CT ,( ri )(5 a s (l)w+ m , CT ,(l) I v± (2)> (75) 

i'm'cr' 

that are needed for evaluating the spherical matrix elements (34) of the Green 
function representing the flow scattered from the walls. 

Recalling notation (2) and definition (16) we note that the basis fields vr mo .(r 2 ) 
are singular at r = R 2 . Accordingly, the transformation formula (74) is valid 
for ±(R 2 — Ri) • e 2 > 0. We also note that the integral defining the matrix 
element on the right side of (74) is not absolutely convergent for / + o < 2, 
because of the slow convergence of the field (16) at infinity; the principal- value 
interpretation of the integral is employed in this case. 

Using the Cartesian displacement relations (66), the matrix elements on the 
right side of equations (74) and (75) can be factorized into the products of the 
displacement matrices (68) and the position-independent matrices T^s (k, lm) 
and Tg(f (Zm, k), 

(5o C (l)w^(l) I vj7 mV (2)> = [S^(R 12 ,k) .T± s -(k,ZW)l , (76) 



(5 a S (l)w+ mCT (l) I v^(2)> 



T+^(/m,k') 



S( . \un. a 1 •Sc ± (Rl2,k') . 

J crcr' 

Equation (68) indicates that the matrix element (76) is nonsingular in the 
limit R12 — > 0, even though the integrand in the scalar product is singular 
at r = for Ri = R 2 . (In contrast, relation (77) does not involve any singu- 
lar integrals.) In the limit Ri 2 — > 0, equations (72) and (74)-(77) yield the 
transformation relations 



(77) 



lmcr\ 



= j dk' ]T v± CT ,(r)T± s -(k', lm; a' | a), ±z < 0, (78) 
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I'm'a' 



v? mV (r)7^(ZW,k;</|<7). (79) 



The matrices T^g (k, lm) and Tg^ (Zm, k) have several important symmetries. 
First, we recall that the Cartesian basis sets (50) and (51) are related to each 
other via the reflection with respect to the plane z = 0. The corresponding 
symmetries of the transformation matrices are 

T cs (k,lm;a\ a') = (-1)'+"^'T+ S (k, lm; 2 - a | a'), (80) 



T+ c -(Zm,k ;0 - | a') = (— l)' +m+cr Tg ( t (/m, k; o" | 2 — a'). (81) 

Another important symmetry relation is associated with the representations 
(21) and (59) of the Oseen tensor in the spherical and Cartesian bases. The 
relation is obtained by applying the Oseen integral operator with the kernel 
in the respective forms (216) and (59) to the fields wj^., which yields 

/ T ( ri -r;)5 V 2 )w^' 2 )dr'= £ v+ mV ,( ri )(v,7 m , CT ,(l) | 5 c (2)w^(2)) 

J I'm'a' 

(82) 

and 

/ T ( ri - r;)5 c (r' 2 )w^(r' 2 ) dr' = J dk'£ v^( ri )(v^(l) | <5 c (2)w± (2)>. 

a' 

(83) 

By comparing the above expressions in the limit R i2 — > we find 

v£(r) = E v+ w (r)T±-*(k, lm; a \ a'), (84) 

I'm'a' 

where equations (72) and (76) and the orthogonality condition (63) were ap- 
plied. Since the expansion of the flow fields v^ cr (r) into v / t m / cr /(r) is unique, 
equations (79) and (84) imply the symmetry 

T£f(Jm,k)= [T± s -(k,Hl f - (85) 



The functional dependence of the matrices Tg^ and T^g" on the wave vector k 
can also be derived using symmetry considerations. Specifically, one can show 
that 



T s + c ± (/m, k; a \ a') ~ k~ 1/2 k l+a ~ l e im ^, (86a) 
T^g" (k, lm; a\a')~ Ar 1 / 2 /^' -1 e" 1 "^, (86b) 
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where k = (k, ip) is the representation of the wave vector in the polar coordi- 
nates. The angular form of relations (86) stems from the requirement in the 
definition of the basis fields (16) and (17) that the coefficients V z ^ lcr (^,0) are 
combinations of spherical harmonics of order m. The dependence on the am- 
plitude of the wave vector k results from the invariance of the transformation 
relations (78) and (79) with respect to the scale change 

r^ar, k -> a _1 k, (87) 

where a is a real parameter (cf., expressions (16), (17) and (50), (51) for 
the spherical and Cartesian basis fields). Using equations (86), the matrices 
TgQ (lm, k) and T^s (k, lm) can be represented in the factorized form 

T+± (lm, k) = (-i) m {2nky 1 / 2 e" im ^ K(k, I) • t+f (Im), (88a) 

T^g(k,H = i m (27rA;)- 1/2 e im ^t± s -(/m)-K(A;,/), (88b) 
where K(k, I) is a diagonal matrix with the elements 

K(k,l;a\a') = 5 aa/ k l+a -\ (89) 



A further simplification of the structure of the transformation matrices Tg^ 
and TqcT results from the curl relations (18)-(20) and (56)-(58). By taking 
curl of both sides of equations (78) and (79), applying the symmetries (80), 
(81), and (85), and using the factorization formulas (88), one can show that 
the matrices Tg^(lm) and T^g (lm) have the following triangular structure 



T++ 
1 sc 



a b c 
2a 2b 
4a 



1 SC 



c b a 
-2b -2a 
4a 



(90a) 



1 cs — 



-1) 



l+m 



c -2b 4a 
b -2a 
a 



a 
6 2a 
c 2b 4a 



(90b) 



and involve only three independent coefficients. As shown in Appendix C, the 
expressions for the coefficients a, 6, c are 
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a = 



[4(/-m)!(/ + m)!(2/ + l)]" 1 / 2 , 



(91a) 



b = 2am/l, 



(91b) 



c — 



a 



l(2l 2 - 21 - 1) - 2m 2 
1(21 - 1) 



(1-2) 



(91c) 



Relations (68)-(73) and (88)-(91) represent the key results of the analysis pre- 
sented so far. In Section 5 we apply these results to express the spherical ma- 
trix elements of the free-space Green operator (33) in terms of two-dimensional 
Fourier integrals (which can be explicitly evaluated in this case). The Fourier 
representations of matrix elements (34) for a system bounded by a single pla- 
nar wall and by two parallel planar walls are derived in §6 and §7, respectively. 
These results enable efficient numerical evaluation of the multiparticle friction 
matrix in wall-bounded suspensions. Description of our algorithm is given in 
Section 8, and examples of numerical results are shown in Section 9. 



5 Fourier representation of the spherical displacement matrix 

In this section we apply the Cartesian displacement formulas (70) and trans- 
formation matrices (88)-(90) to express the spherical displacement matrix 
Sg _ in terms of two-dimensional Fourier integrals. Such a representation can 
be utilized in developing numerical algorithms for evaluating hydrodynamic 
interactions in doubly periodic systems. Moreover, the analysis allows us to 
introduce some techniques that will be used in the discussion of the flow in 
wall-bounded suspensions (cf. Sections 6 and 7). 

We recall that the displacement matrix Sg _ and the corresponding spherical 
matrix elements (33) of the Oseen operator are equivalent, as indicated by the 
formula (43). To make the notation in this and the following sections parallel, 
we express our results in terms of the matrix elements G° (/m | I'm'). 

By inserting the expansion (78) into (41) and using equation (77) we find 



where the plus sign applies for R^- • e 2 < and the minus sign for R^- ■ e z > 0. 
We note that the Lorentz symmetry (45) of the matrix elements G° (/m | I'm') 
is explicit in equation (92) due to the symmetry relations (73) and (85) for the 
component matrices. The angular integration in relation (92) can be explicitly 
performed with a help of the factorization formulas (68) and (88) and the 




(Im | I'm') = rj 



-i 



/ 



dkT+ c ± (/m, k) • S± ± (R i , k) • T± s "(k, I'm'), (92) 
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equation 

f ^ e Kkp cos ^-mip) d ^ = 2ni m J m (kp), (93) 

Jo 

where J m (x) is the Bessel function of the order m. The resulting expression 
has the form 



G%(lmo | I'm! a') = rj' 1 (-l) m ~ m e^ m ~ m) ^ 

x / gtikZ^; Ima \ l'm'a')k l+l + ° + ° - 2 J„'(^) dk, (94) 

J 



where (pij, 4>ij, Zij) is the representation of the vector Ry in the cylindrical 
coordinates, and 

g^(kz;lm | I'm') = t+f (Zm) • S^faz) ■ (I'm'). (95) 



Relations (70) and (95) indicate that the integrand in equation (94) is a com- 
bination of the Bessel function, powers of k, and the exponential e~ fc ' z ^L The 
integrals of this form can be evaluated using the following identity 

k l J m (kp)e- kz dk = (-l) m (Z-m)!r- (m) PJ n (r- 1 z), z > 0, (96) 

o 

where PJ n (x) is the associated Legendre polynomial, and 

r = (/ + z2) i/ 2 _ (Q7) 

We have verified that equations (94)-(96) yield expressions that are equivalent 
to the displacement theorems for the spherical basis of Stokes flows derived 



by Felderhof and Jones [411 ] . In development of the multipolar-expansion al- 
gorithms to evaluate hydrodynamic interactions in doubly-periodic systems, 
a direct application of the Fourier representation (92) may be useful. 



6 Single-wall Green operator 



Similar techniques can be used to evaluate the matrix elements of the Green 
operator (34) in a system bounded by a single planar wall. We assume that 
the wall is in the plane 

z = Z w (98) 
and the suspension occupies either the halfspace z > Z w (denoted by + ) 
or z < Z w (denoted by The spherical matrix elements of the Green 

operator (34) for this system are obtained by combining the transformation 
relations (88)-(90) between the spherical and Cartesian basis sets with the 
Cartesian representation of the flow reflected from the wall. The reflected flow 
is discussed in the following subsection. 
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6.1 Single-wall reflection matrix 



The velocity field in the halfspace f^, occupied by the fluid, can be uniquely 
decomposed into the incoming and reflected flows 



v r 



,out 



where 



(99) 

(100a) 
(100b) 

(101) 

denotes the position of the point r relative to the wall, where R w = (X w , Y w , Z W/ 
has arbitrary lateral coordinates X w and Y w . 



= I dk^C(ka)v±(r w ), 

J a 

vr(r) = / dk^C 4 Nv£(rw). 



In the above relations 



According to the definitions (50) and (51), the decay of the basis flow fields 
v kcr( r w) for k — > oo is faster in the halfspace Q^ 1 than it is on the wall surface 
(98). Thus, assuming that the integral (100b) converges on the surface (98), 
the scattered flow field v° ut (r) is nonsingular in the whole region fi ± occupied 
by the fluid. Likewise, the convergence of the integral (100a) on the wall surface 
implies that the incoming flow field v™(r) is nonsingular in the complementary 
region Q T . 

By analogy with the relations (39) and (40) for a flow field scattered by a 
particle, we introduce the single- wall scattering matrix Z w , defined by the 
equation 

cr(k) = -Z w -C(k), (102) 

where c° ut (k) and c™(k) denote the arrays of expansion coefficients in equa- 
tions (100). For an immobile rigid wall, the velocity field (99) vanishes at 
z = Z w . By inspection of expressions (50) and (51) we find that 



1 
1 
1 



(103) 



in this case. For planar interfaces with other boundary conditions (e.g., a 
surfactant-covered fluid-fluid interface discussed in 42||) the scattering matrix 
is different from identity, and it may depend on the magnitude of the wave 
vector k. Explicit expressions for scattering matrices for such systems will be 
derived in forthcoming publications. 
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6.2 Matrix elements of Green operator 



In order to evaluate matrix elements (34) of the single-wall Green operator we 
consider the flow field produced by the induced-force distribution (1) centered 
at the position of particle j. By comparing the decompositions (4) and (99) 
and using relation (78) we find 

= /T (r-r').F,(r')dr' (i 4) 

and 

vrW= /T'(r,r')-F,(r')dr'. (l 05 ) 

We note that according to equations (11) and (104) the flow incident to the 
wall equals to the flow scattered by the particle 

v»=vf(r). (106) 

By projecting equation (105) onto the reciprocal spherical basis w^ CT centered 
at the position of particle i and using the multipolar expansion (29) we get 

<<S(0wjL(0 I vD = E G'^lma | I'm'a') (I'm' a'), (107) 

I'm'a' 

where the definition (34) was applied. The matrix element of the reflected flow 
at the left side of the above equation is evaluated with the help of the reflection 
relation (102). Accordingly, the expansion coefficients of the incoming flow 

C(k<r) = (5^(w)wt(w) | O (108) 

(where the index w in the bra (5q (w)w^ a (w) \ indicates the dependence on the 
variable (101)) are determined using expansion (40) for the incoming velocity 
field (106) and the relation (76) for the matrix elements relating the spherical 
and reciprocal basis fields. Collecting these formulas yields 

C(k) = r)- 1 J2 S^(Rwi, k) • T± s -(k, I'm') • fjil'm'), (109) 

I'm' 

where Rj w = Rj — R w and R WJ - = R w — TLj. The above relation is combined 
with the expansion (100b) and the scattering formula (102); the resulting 
expression for v° ut is inserted into (107). The matrix elements between the 
spherical and Cartesian basis fields are then evaluated using (77). By compar- 
ing the result of this calculation to the expression (107) we find 

G'ijilm | I'm') = -rj- 1 ! dk T+f(Zm, k) • S^(R lw , k) • Z w • S± ± (R wj , k) • T± s (k, 

(110) 
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A physical interpretation of the above relation is straightforward: the spherical 
components of the flow produced by the particle at point j are transformed 
into the Cartesian basis by the matrix T^g"; the Cartesian components are 
propagated by the matrix Sc ± (R w j) to the wall, where they are scattered 
(as represented by matrix Z w ); the reflected field is propagated to the point 
i by the matrix SQ T (Rj w ); and finally the flow is transformed back into the 
spherical basis by the matrix Tgjf . We note that, similar to relation (92), the 
Lorentz symmetry (46) of the matrix elements (110) is explicit due to the 
symmetry relations (73) and (85) of the component matrices. 

Similar to the angular integral in equation (92), the angular integration in the 
Fourier representation (110) of the matrix G'jj can be explicitly performed 
with the help of expressions (68), (88), and (93). The resulting expression for 
the matrix elements of the one-wall Green operator is 

G'^ilmo | I'm' a') = r ] - 1 (-l) m '- m e i( - m '- m ^» 

x g^{kZ iW) kZ wr M<y\l'm'a')k l+l+ ^- 2 l m ^ m {kp l3 )^ (111) 

J 

where 

g±(kZ iw , kZ wj] Im | I'm') = -T&f (Zm) • S^(kZ iw ) • Z w • S^(kZ WJ ) • f ^(ZW). 

(112) 

We recall that the upper signs in the above equations correspond to the fluid 
occupying the upper half-space Q + , and the lower signs to the fluid in the lower 
halfspace Q~. Taking this into account, we find that the exponential factors 
resulting from the Cartesian displacement matrices (70) in the product on the 
right side of equation (112) can be combined into a single exponential factor 

g±(kZ iw , kZ WJ ; Im | I'm') ~ e- fcA ^ (113) 

where 

Ay = \Z lw \ + \Z wj \ (114) 

is the vertical offset between the target point i at the position Rj and the 
image of the source point j at 



Rj = Rj - 2(Zj - Z w )e z . (115) 

It follows that the integrand in equation (111) is the combination of the fac- 
tor (113), the Bessel function, and powers of k. Thus, relations (96) and (97) 
imply that the elements of the matrix can be expressed in terms of the 
flow produced by an image singularity at r = RA. We note that such a form 
is required by the Lorentz reflection relation 43||. Explicit expressions for the 



image force multipole system corresponding to an arbitrary source force mul- 
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tipole have recently been derived by Cichocki and Jones [20(; we have verified 
that the integral (111) yields results equivalent to their expressions. 



The main application of the Fourier representation of the single-wall Green 
operator G'jj is in the subtraction technique that is implemented in our al- 
gorithm for evaluating the multiparticle friction matrix in a two-wall system. 
In this application (in more detail described in Section 8.2) expressions (111) 
and (112) are used in conjunction with the results of Ref. |20j to accelerate 
the convergence of the Fourier integrals for the two-wall Green operator. 



7 Two-wall Green operator 

In this section we generalize the analysis of Section 6 to a system of particles 
confined between two parallel planar walls. We assume that the walls are in 
the planes 

z = Z L , z = Z v , (116) 

where 

Z L < Z v . (117) 
The suspension occupies the region Zl < z < Zjj. The positions of the walls 
are indicated by vectors Rl = (Xl, 1l, Zi) and Ru = (X\j,Y\j, Zjj), where 
the lateral coordinates Xl and YL and X\j and Y\j are chosen arbitrarily. 

The flow produced in this system by the induced-force distribution (1) centered 
at the position of particle j is a superposition of three components 

v(r) = v° ut (r) + v° ut (r) + v° ut (r). (118) 

Here v° ut (r) is the velocity field (11) produced by force distribution Fj, and 
v° ut (r) is the flow reflected by wall a = L, U. By definition, the flow component 
v£ ut (r) is nonsingular in the region z > Z^ and vanishes for z — > oo, and the 
flow component v™*(r) is nonsingular in the region z < Z\j and vanishes for 
z — > — oo. Accordingly, the expansions of the flow fields v£ ut and v™' in the 
Cartesian basis sets (50) and (51) have the form 

v£*(r) = / dk$>r(ka)v k(7 (r L ), (119a) 

v° ut (r) = / dk]Tcr(k<7)v+ ( ru)j (ugh) 

where tl = r — Rl and ru = r — Ru- Expressions (119) are consistent with 
the expansion (100b). 

The three components (118) of the velocity field produced in the space be- 
tween the walls by the force distribution Fj can be used to construct the flow 
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components (a = L, U) incoming to wall a. Using expressions (119) and 
the definition (100a) of the incoming flow we find 



vL n (r)=v° ut (r) + v° ut (r), 



(120a) 



v^rHvrM + vHr). 



(120b) 



Relation (11) and the respective decompositions (4) and (118) of the Green 
function T(r, r') and the flow field v(r) imply that 



vT t (r)+vS rt (r) = /r^rO-F^rOdr'. 



121) 



Projecting the above equation onto the reciprocal spherical basis w^ CT yields 
(<^)w+ mCT « I O + ^Ww^W | O = £ ^.(W | l'm'a')J)(l'mW), 

I'm' a 1 

(122) 

which is analogous to the single-wall result (107). Explicit expressions for the 
matrix elements G'^ilma \ I'm' a') can thus be derived by generalizing the 
analysis presented in Section 6. 

To this end, the representation of the velocity fields v™ in terms of the Carte- 
sian basis fields v^ : (T (r Q ,) aligned with the wall a is obtained by inserting ex- 
pansions (40) and (119) into (120), and applying the transformation formulas 
(66) and (76). Using then the single- wall scattering formula (102) to relate the 
expansion coefficients for the outcoming and incoming flows we get a pair of 
coupled linear equations 

c° ut (k) = -Z w • [S++(R LU , k) • cr^+TT 1 E SJ + (Rl„ k) • T+ s (k, I'm') ■ f(l'm% 

I'm' 

(123a) 



cT (k) = -Z w • [S5"(R UL , k) • cr (lO-HT 1 E S c ~(Ru„ k) • Tcs (k, I'm') • f (I'm')}, 

I'm' 

(123b) 

R a , /3 = R a -R /3 £ = L,U,j. (124) 



where 



In order to to express the solution of the system (123) in a compact manner 
we introduce a matrix notation in the space of six-dimensional vectors of the 
form 

'cHk) 



c(k) = 



cr(k) 



(125) 
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Accordingly, we define the 6x3 and 3x6 transformation matrices 



T cs (k, lm) 



Tcs ( k , lm ) 
T cs (k, lm) 



T sc (lm,k) = 



T+ c (/m,k) T s + +(/m,k) 



the 6x6 Cartesian displacement matrices 



Swj(k) 



S+ + (R Lj ,k) 



S c -(Ru„k) 



(126a) 



(126b) 



(127a) 



S5"(R lL ,k) 

Siw(k) = 

S++(R,u,k) 
and the 6x6 two-wall flow reflection matrix 



Z TW (k) 



S++(R LU ,k) 



S c -(R UL ,k) 



r-l 



(127b) 



(128) 



For simplicity, the dependence on the wall and particle positions was sup- 
pressed on the left side of the above expressions. 

Due to the symmetries (73) and (85) of the 3x3 transformation and displace- 
ment matrices, the corresponding symmetry relations 



Tcs(k,/m) = [7 S c(/m,k)]t, 



(129a) 



S w ,(k) = [S iW (k)] t , 



(129b) 



Z TW (k) = [Z TW (k)]t (129c) 

are satisfied by the matrices (126)-(128). We note that the two- wall scattering 
matrix (128) involves displacement components describing translation of the 
flow field between walls. 
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Using notation introduced above, the solution of the system (123) can be 
written in the form 



c(k) = -rf 1 £ Z TW (k) • S w ,(k) • 7 C s(k, I'm') ■ f,-(ZW). 



(130) 



I'm' 



In order to get an explicit expression for the matrix elements of the Green 
operator G'y, equations (119), (122), and (130) are combined, and the scalar 
products between the Cartesian and spherical basis fields are evaluated with 
the help of relation (77). The resulting expression for the elements of the 
two-wall Green matrix is 

G'ijilm | I'm') = -rT 1 / dkT sc (lm, k) • S lW (k) • Z TW (k) • S w ,(k) • 7 cs (k, I'm'). 

(131) 

The expression is similar in its form (and interpretation) to the corresponding 
relation (110) for a single-wall system. As with the matrix elements (92) and 
(110), the Lorentz symmetry (46) of the elements (131) is explicit due to the 
symmetry relations (129) of the component matrices. 

The dependence of the integrand in equation (131) on the polar angles in 
the Fourier and real spaces is identical to the corresponding dependence in 
equation (110) — the angular integration can thus be performed in a similar 
manner. By analogy with equations (111) and (112) we get 



G'^ilma | I'm' a') = r]' 1 {-l) m ' ~ m jW-™)** 

roc 

x J gTw(k; Ima \ l'm'a')k l+l +a+a ~ 2 J m ^ m (kp i:j ) dk, (132) 

where 

g TW (A;; Im \ I'm') = - t sc (lm) • S m (k) • Z TW (k) • S Wj (k) • t cs (l'm'), (133) 
with the matrices in the product given by 



t cs (lm) = [t sc (lm) 



T cs ( lm ) 



t cs (lm) 



(134) 



S Wj (k) = [S jW (/c)] f = 



S+ + (kZ Lj ) 



Sc~ (kZ Vj ) 



(135) 
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and 



Z TW (k) 



where 



Z-J S+ + (-kH) 
H = Zij — Zj 



(136) 



(137) 



is the wall separation. The above relations, together with equation (70), im- 
ply that gTw(&; l m I I'm!) depends on the z-coordinates of the walls and the 
points i and j, but is independent of the lateral coordinates, consistently with 
the translational invariance of the system. Since the two-wall scattering ma- 
trix (136) is more complex than its one-particle counterpart, the integration 
in equation (132) (unlike (111)) cannot be analytically performed. However, 
numerical integration is straightforward, except when the lateral distance be- 
tween points i and j is large, in which case the oscillatory character of the 
integrand becomes important. 



8 Stokesian-dynamics algorithm for suspension between two walls 



The theoretical results derived in the previous sections enable development 
of efficient numerical algorithms for evaluation of many-body hydrodynamic 
interactions in suspensions of spherical particles confined between two planar 
walls. To our knowledge, such algorithms have not been available so far. In 
what follows, we describe a many-particle Stokesian-dynamics algorithm based 
on our theory. 

In Section 8.1 we summarize expressions that relate the matrix M in the 
force-multipole equation (48) to the resistance matrix in a suspension of many 
spheres. Our transformation formulas relating the spherical and Cartesian ba- 
sis fields are employed in Section 8.2, where a simple numerical-integration 
procedure for evaluating the elements of the matrix M is described. The 
lubrication-subtraction techniques [15. Il7l l2l| used for improving convergence 



with the order of the force multipoles included in the calculation are outlined 
in Section 8.3. Examples of numerical results for the friction matrix of a sin- 
gle particle, a pair of particles, and in many-particle systems are presented in 
Section 9. 
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8.1 Resistance matrix 



We focus on a system of N spheres undergoing translational and rotational 
rigid-body motion (7) with no external flow. The particle dynamics in the 
system is characterized by the resistance matrix 



c 



IJ 



Ctt >tr 
ij Sij 

c rt c 



rt >rr 
ij Sij 



defined by the linear relation 



hJ 



,N, 



(138) 



N 

E 

3=1 



■-tt >tr 



Crt >rr 



(139) 



between the translational and rotational particle velocities TJj and flj and the 
forces and torques Ti and Xj applied to the particles. The dot in equation 
(139) denotes the matrix multiplication and contraction of the Cartesian ten- 
sorial components of the resistance matrix. A detailed discussion of a more 
general resistance problem, which involves an external linear flow and the 



stresslet induced on the particles, is presented in Ref. |21 



The resistance relation (139) can be linked to the induced-force equation 
(30) by expressing the applied forces and torques 7~i and T*j in terms of 
the induced-force distributions (1), 



Fi = J P f (r) dr, Ti = J r^Fifr) dr. 



(140) 



Representing the above quantities in terms of the induced-force multipoles 
(29) yields 

Fi = J2 X (t I Ima) f (Ima), T, = X(r | Ima) f (Ima), (141) 

Ima Ima 



where 



X(t | Ima) = fc^XV) 



142a) 



X(r | Ima) = SnS^X^m). (142b) 

Explicit expressions for the vectors X t (m) and X r (m) are listed in Appendix 
B. The coefficients Cj in the corresponding expansion (36) for the flow field (7) 
can be represented in the form 

d(lma) = X(lma | t) • U» + X(lma \ r) • n t , (143) 
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where v ext (r) = is assumed. As shown below, the projection matrices X in 
equations (141) and (143) obey the identity 



X(lma | A) = X*(A \ Ima), A = t,r. (144) 



In order to determine the resistance matrix £y from the above expressions, 
the force-multipole equation (48) is solved, which yields the linear relation 

N 

f^lm) = Y,Y, F ij( lm I l ' m ') • Cj(/W), (145) 

j=l I'm' 

where F = M _1 is the generalized friction matrix. By inserting into equation 
(145) the projection formulas (141) and (143) we get 

Cfj^ = J2 £ X ( A I lma)Fij{lma \ l'm'a')X(l'm'a' \ B), (146) 

Ima I'm' a' 

where A, B = t, r. 

With our normalization of the spherical basis fields (16) and (17) (as defined 
in Section 3.1) the symmetry relation (144) is a direct consequence of the 
Lorentz symmetry of the generalized friction matrix 

Fijilma | I'm' a') = Fji(l'm'a' | Ima), (147) 



which follows from equations (45)-(47). Relation (144) is obtained by inserting 
equation (147) into (146) and using the Lorentz symmetry of the resistance 
matrix [14J 



where the dagger denotes the transposition of a tensor. 



8.2 Evaluation of matrix M 



The evaluation of the resistance matrix C,- from expression (146) requires 
solving the set of linear algebraic equations for the induced-force multipoles 
(48) to obtain the generalized friction coefficients Fijilma \ I'm 1 a'). The matrix 
M in the equation (48) is given as the sum of three terms (49). The first two 
terms, i.e., the single-particle scattering matrix Zj and the matrix G°- repre- 
senting the flow in infinite space, are known explicitly 0,E3|. The remaining 
term — the two-wall contribution G'ij — is evaluated numerically, using rela- 
tions (132)-(136) along with our expressions for the Cartesian displacement 
matrices (70), the transformation matrices (90), and the single- wall scattering 
matrix (103). 
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k k 

Fig. 1. Integrands * and 5^, denned by equations (149) and (150), versus magni- 
tude of the wave vector k. Separation between walls H = 1; distance of the source 
and target points from lower wall hi = hi = 0.1. 

As already mentioned at the end of Section 7, numerical evaluation of the inte- 
gral (132) is straightforward for sufficiently small lateral separations between 
particles i and j. For large interparticle separations p^-, however, the inte- 
gration is more difficult because of the oscillatory behavior of the integrands 

= 9Tw(k; Ima | lW)k l+l ' +a+a '- 2 J m ,- rn (kp ij ), (149) 

due to the presence of the factor J(kpij). This behavior is illustrated in the 
left panel of figure 1 for a configuration in which both points % and j are close 
to one of the walls. 



To avoid numerical integration of a highly oscillatory function, the integrand 
(149) is decomposed 

^(k) =^ L (k) + ^ v (k) + 5^(k) (150) 

into the superposition of the single-wall contributions \&l and \I/tj, an d the 
wall-interaction part 5^/. According to equations (111) and (112), the single- 
wall integrands are 

V a (k) = g±(kZ ia , kZ af , Ima | I'm! o')k l + l '+°+°'- 2 l m ,_ m {kp l0 ), (151) 

where a = L, U. Relation (113) implies that for large values of k the amplitude 
of these integrands decays as 

* a (k)~e- kA( Z\ (152) 

where is the vertical offset (114) between the point i and the reflection 
of point j in the wall a. Therefore, the decay is slow if both points i and j are 
close to the wall, consistent with the results in the left panel of figure 1. 

In contrast, the decay of the wall-interaction part 5^/(k) of integrand (150) is 
determined by the wall separation H. As shown in Appendix D, the large-/c 
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asymptotic behavior of this function is 

6V(k) ~e" fcA ^, 



(153) 



where 

Ay = 2H- \Zij\ > H. (154) 

The lengthscale A^- equals the vertical offset |Zj — Z"\ between the target 
point i and the closer of the two second-order images of the source point j. 
The images are at the positions 

RJ = R,- ± 2iZe*, (155) 

where the plus sign corresponds to reflecting the source point first in the lower 
wall and then in the upper wall; the minus sign corresponds to the opposite 
order of reflections. 

A typical form of the wall-interaction contribution S^f(k) is presented in the 
right panel of figure 1. Unlike the results for the integrand 5^f(k) in 

this example is negligible already after several oscillations. Thus, the function 
5ty(k) is easy to integrate numerically. 

In our algorithm, the contribution SG'^ilma \ I'm' a') to the matrix elements 
(132), associated with the component 5^(k) of the integrand, is evaluated 
by numerical integration using the Simpson rule. The slowly convergent one- 
wall contributions (151) are calculated analytically, using the explicit image- 
representation expressions (cf., the discussion in Section 6.2). 

Our numerical tests indicate that this procedure yields accurate results for 
Pij/H < 20. The procedure can be further improved, either by subtracting 
several terms associated with higher-order wall reflections of the source mul- 
tipole [221, or by deriving asymptotic formulas for the integrals (132). 

8.3 Convergence with multipolar order 

Our numerical procedure for evaluating the friction matrix involves truncation 
of the linear system (48) by neglecting the induced force multipoles fj(Zm) of 
the order / > Z max . This multipolar approximation converges very slowly with 
the truncation order Z max if any two particles are close together or a particle is 
close to a wall. Such a behavior stems from a rapid variation of the flow field 
in the near-contact lubrication regions. The mechanism is well known and has 
been observed for particles in infinite space and in the presence of a single 
wall. To overcome this difficulty we employ a standard method, originally 
introduced by 0], according to which the lubrication forces are included in 
the friction matrix using a superposition approximation. We follow closely the 
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Fie. 2. Relative error 5 a of trie lateral (^||) and vertical components of the trans- 
lational friction matrix (160) for a single sphere between two parallel walls, versus 
truncation order l max in the multipolar approximation (159). Left panels correspond 
to center and right panels to off-center particle position (as indicated). Values of 
dimensionless gap (162) between the particle and the closer wall are e w = 0.02 (open 
triangles); e w = 0.1 (circles); e w = 1.0 (solid triangles). 



implementation of the method described by |21| in their study of a single wall 
problem. Accordingly, the resistance matrix (138) is represented in the form 



t = /* SU P' 2 _|_ A SU P> W _i_ A/- 



where 



and 



c 



sup, 2 
ij 



N 

'•3 E 

k = 1 
k^i 



h 3 £ C {0) ^^(^k) + (l-5 lJ )C [J M^J) 



(o), 



c 



sup.w 



Said 



(156) 



(157) 



(158) 



Here Cmm( mn ) * s the se h~ an d Cmn( mn ) the mutual- resistance matrix for an 
isolated pair of particles m and n in the unbounded space, and Cm( m ) * s the 
single-particle resistance matrix for a sphere in the subspace bounded by the 
wall a = L,U. The superposition contributions (157) and (158) in equation 
(156) can be calculated with high accuracy, using methods discussed below. 
The convergence with the truncation order Z max of the multipolar approxima- 
tion 



■.sup, 2 



+ c p,w + t A c 



I J Umax 



(159) 
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is much faster than the convergence of the multipolar approximation [C^luax 
itself, where [5]/ max denotes the quantity B evaluated using the multipolar 
expansion truncated at / = / max . Therefore the evaluation procedure based on 
equation (159) yields accurate results for the friction matrix at a substantially 
reduced numerical cost. 

In our implementation, the two-particle components (^ m (mn) and (^(mn) 
of the friction matrix in the superposition formula (157) are evaluated using 
a combination of the lubrication resistance expressions 3 and the series 



expansion in inverse powers of interparticle separation |3Jj. Similarly, the one- 
particle friction matrix Cm( m ) i n the superposition formula (158) is evaluated 
using a combination of the lubrication resistance expression and the power 
series in the inverse distance of the particle from the wall |2o| . 

Our numerical results indicate that for large and moderate wall separations H 
(compared to the particle diameter) the multipolar approximation (159) con- 
verges rapidly with the truncation order Z max . For configurations with H ^ 2a 
the convergence is less satisfactory, particularly for the transverse components 
of the friction matrix. This behavior is illustrated in figure 2, where the relative 
error for the lateral and vertical components 

di = Cn xx = err, a = ct\ zz (m 

of the one-particle translational friction matrix is shown for different particle 
configurations and distances between the walls. The results are given for the 
center and an off-center position of the particle in the space between the walls, 

h = \H, h = \H, (161a, b) 

where h is the distance of the particle from the lower wall. In the case of 
the center particle position (161a), the multipole-truncation error exhibits 
decaying oscillations. For small values of the gap 

e w = h/a-l (162) 

between the particle surface and the closer wall a typical error is in the range of 
several percent. The results corresponding to truncations at even orders of Z max 
are more accurate than the results corresponding to truncations at odd orders. 
The multipolar-truncation error for the off-center particle position (1616) is 
much smaller than the error for the center configuration with the same values 
of the particle-wall distance h. 

A similar dependence of the truncation error on Z max was observed for many- 
particle friction matrix: the error is small, except when the wall separation is 
comparable to the particle diameter. This behavior suggests that the relatively 
large error for such tight configurations results from an interaction between the 
lubrication layers — this effect is not accounted for in the superposition terms 
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Fig. 3. Normalized translational resistance coefficients (163) for a pair of particles 
on the x axis in the mid-plane between the walls, versus the interparticle distance 
P12 normalized by particle diameter d. Wall separation H/d = 1.1 (solid lines); 2.0 
(dashed lines). Heavy lines represent the exact results and thin lines the superposi- 
tion approximation (174). 

in equation (159). The problem, however, requires further investigations in 
order to develop better approximation methods. 



9 Numerical results 



In this section we give some examples of numerical results for the hydrody- 
namic friction matrix in systems of spherical particles confined between two 
parallel planar walls. The calculations for a single particle and for particle 
pairs, depicted in Figs. 3 and 4, were performed using the multipolar approx- 



34 




Fig. 4. Normalized one-particle translational resistance coefficients (167) for a par- 
ticle at the center position in the space between the walls, versus wall separation H 
normalized by particle diameter d. Horizontal component C, = Cii (solid lines); verti- 
cal component C = C± (dashed lines). Right panel represents the relative accuracy 
C s /C of the superposition approximation (174). 



imation (159) with the truncation at the order Z max = 12. The multi-particle 
calculations, depicted in Figs. 5 and 6, were obtained using Z max = 8. These 
truncations are sufficient to obtain results with the accuracy better than the 
resolution of the plots. A more extensive set of numerical results is presented 
in a separate publication [44 . 



9. 1 Two-particle friction matrix 



Figure 3 illustrates the behavior of the translational components of the two- 
particle resistance matrix, normalized by the Stokes friction coefficient Co — 
6m]a, 

C,* C7Co. «,i = l,2, (163) 

where a,/3 — x,y,z. The particle pair is in the center plane of the space 
between the walls 

ht = h 2 = \H, (164) 

where hi is the distance of particle i from the lower wall. The relative particle 
displacement is along the x direction 

P12 = P12&X, (165) 

and the results are plotted versus the interparticle distance pi 2 . Only the 
diagonal Cartesian components Cn* an d C12 °f the self- and mutual- resistance 
matrices are shown, because S = for a 7^ (3, due to symmetry. 
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4.0 





Fig. 5. Translational resistance coefficients per particle (173) of rigid linear arrays 
of touching spheres on a line parallel to axis x at the center plain between the 
walls, scaled by corresponding one-particle values (167), versus the wall separation 
H normalized by particle diameter d. Number of spheres N = 2 (solid line); 5 
(dashed); 10 (dot-dashed); 20 (dotted). 

9.1.1 Near-contact and intermediate behavior 



According to the results shown in the left panels of Fig. 3, the self-components 
of the two-particle resistance matrix are only weakly affected by the presence 
of the second particle, except for sufficiently small gaps between the particle 
surfaces 

e = p 12 /d-l, (166) 

where d = 2a is the particle diameter. The effect of the interparticle interac- 
tions is most pronounced for the longitudinal component £ff , because of the 
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Fig. 6. Same as Fig. (5), but for the relative accuracy Cc aS /Cc Q °f superposition 
approximation (174). 

strong 0(e _1 ) lubrication resistance for particles in relative motion along the 
line connecting their centers. For the motion in the transverse directions y 
and z, a significant interparticle-interaction effect is seen only for very small 
interparticle gaps, because the transverse interparticle lubrication resistance 
has a much weaker logarithmic singularity O(loge) than the longitudinal one. 

The results in the right panels of Fig. 3 indicate that for small interparticle 
distances all three components £f 2 a , a = x,y,z, of the mutual friction matrix 
are negative. The negative sign indicates that the hydrodynamic force J-^ = 
—J-\ produced on particle (1) by the motion of particle (2) points in the 
same direction as the particle velocity. When the distance between particles 
is increased, the transverse components Cff and £f| change sign, which results 
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from the backflow associated with the flow field scattered from the walls. In 
contrast, the longitudinal component £ff remains negative. We note that the 
backflow effect does not occur in the unbounded space. 



9.1.2 Far- field behavior 

At large interparticle separations pxz/d ^> 1, the mutual components °f 
the two-particle resistance matrix vanish, and the self-components Cn* tend 
to the corresponding one-particle values 

C|| = Cll/Co, bar( ± = a/Co (167) 

(Cll for fff and Qx, and bar(^ for Qx)- The lateral and transverse one-particle 
friction coefficients 0| and bar(^ are shown in the left panel of Fig. 4 versus 
the wall separation H for the center particle position h = \H. Our present 
one-particle results agree with those of Jones [3l| and with our earlier solution 
obtained using the image- singularity technique [271 ]. 

The asymptotic approach of the two-particle friction matrix to the limiting 
values at large px2 can be determined from the far-field behavior of the flow 
field v as produced by a particle moving in the horizontal direction e a (a = 
x, y). By expanding the flow v as in the small parameter H/ p, where p 3> H is 
the horizontal distance from the moving particle, we find that 

v as = \rf x z{E - z)Vp as . (168) 

Here z = is the position of the lower wall, and 

P as ~^, a = x,y, (169) 
p 2 

is the pressure field that depends only on the lateral position p = xe x + ye y . 
The above equation indicates that the far- field velocity (168) decays as 

v as ~ p~ 2 (170) 

for p — > oo. One can also show that the flow field v as produced by a particle 
moving in the z direction decays exponentially. The above results are consis- 
tent with the asymptotic expression for the flow field produced in the space 
between the walls by a Stokeslet pointing in the horizontal direction |45|. (A 
general analysis of the far- field flow will be presented elsewhere.) 

The result (170) implies that the asymptotic far-field behavior of the lateral 
components of the friction matrix is 

CiT = Cn + o( P r 2 4 ), (i7i) 
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CiT = 0( Pu 2 ), (172) 

where a = x,y (the 0(pj~ 2 4 ) contribution in Eq. (171) corresponds to the flow 
field (170) scattered back to the original particle). In contrast to the results 
(171) and (172), the limits £fi = bar(± and Cif = for the vertical components 
of the friction matrix are approached exponentially on the lengthscale H. The 
numerical results shown in Fig. 3 agree with the above analysis. In particular, 
the signs of the longitudinal and transverse friction coefficients and Cif f° r 
Pu/d^ 1 are opposite, consistent with expressions (168) and (169). 

9.2 Linear arrays of spheres 

In order to illustrate the role of the far-field flow for in wall-bounded systems, 
we present, in Fig. 5, the resistance function of rigid linear arrays of N touching 
spheres. The spheres are placed in the mid-plane between the walls on a line 
pointing in the x direction. The figure shows the diagonal components of the 
translational resistance matrix of the array treated as a single rigid body 

N 

Cc° = (N(o)- 1 £ a = x,y,z. (173) 

The normalization of the resistance matrix (173) corresponds to the hydrody- 
namic friction evaluated per one sphere. The results shown in the Fig. 5 are 
further rescaled by the corresponding one particle results (167), and they are 
plotted versus the normalized wall separation H/d for several values of the 
chain length N. 

The results in Fig. 5 indicate that for large separations between the walls 
(compared to the chain length) all three components of the resistance matrix 
(q decrease monotonically with N. Consistent with the behavior of elon- 
gated particles in unbounded space j4(| |4?| , we find that Cc" ~ V 1°§ N 
and (q 1 ~ (q ~ 2(q x for 1 < iV < H/2a. In contrast, for moderate and 
small values of the wall separation H the behavior of each component Cc° °f 
the resistance matrix (173) is qualitatively different. The longitudinal compo- 
nent (q x decreases monotonically with N, while the other two components CcT 
and increase with N due to backflow associated with the presence of the 
walls. This effect is particularly pronounced for the transverse component (q 1 , 
where the increase is by a factor greater than three for iV = 20 in the regime 
H/dm 1.2. 

The qualitatively different behavior of the transverse resistance coefficients 
and C^q is associated with the opposite directions of the asymptotic flow 
field (168) on the horizontal lines parallel and perpendicular to the velocity 
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of a particle. According to relations (168) and (169), the flow field v as in 
front and behind the moving particle points in the direction of the particle 
velocity. This results in a cooperative effect leading to a reduced resistance per 
particle for the longitudinal motion of an array. The direction of the flow on the 
perpendicular line is opposite, which produced a cumulative effect leading to 
a large increase of the resistance coefficient (<j '. This effect is further discussed 
in Ref . 44 1 . 



9.3 Superposition approximation 



To illustrate the effect of the hydrodynamic interactions between walls on 
particle dynamics in wall-bounded systems, the results of our accurate numer- 
ical calculations are compared to the single-wall superposition approximation 

ci c!) • c)j-c:; : . (174) 

In the above equation, (£^) represents the friction matrix for a system of iV 

particles in the presence of the lower (upper) wall, and Cif is the correspond- 
ing friction matrix in the absence of the walls. We emphasize that, unlike the 
superposition terms in equations (156)- (158), all quantities on the right side 
of equation (174) represent the full iV-particle friction matrices — the superpo- 
sition refers only to the wall contributions. The subtraction of the free-space 
term Qj assures that the matrix has a correct limit if the distance of the 
particles to one of the walls (or both walls) tends to infinity. 

In Fig. 3, the translational friction coefficients 

fact S /-tt aa S / /• ( T7K\ 

Uj — Sij /SO, U'OJ 

evaluated in the superposition approximation, are plotted along with the accu- 
rate results for the two particle system. The right panel of Fig. 4 represents the 
ratio CiT/CiT S for 

a single particle. The results indicate that the superposi- 
tion approximation is quite accurate for the single-particle friction coefficients 
and the self-components of the two-particle friction matrix — the maximal er- 
ror for these quantities is about 18% (also see the more extensive one-particle 
calculations reported in Ref. |3l|). 

The superposition approximation is much less accurate for the mutual com- 
ponents °f the two-particle friction matrix, as shown in the right panels of 
Fig. 3. The accuracy of approximation (174) is especially poor for the trans- 
verse component £{* 2 a for small values of the wall separation H. Moreover, the 
approximation yields an incorrect 0(p~ l ) asymptotic behavior of the mutual 
friction coefficients for large p. 
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The failure of the superposition approximation is particularly pronounced for 
the transverse component of the friction coefficient (173) for rigid chains of 
spheres. As shown in the second panel of Fig. 6, relation (174) grossly under- 
estimates the coefficient for long chains, especially for small and moderate 
values of the normalized wall separation H/d. The superposition approxima- 
tion is insufficient in this regime, because it does not accurately reproduce the 
far-field interparticle interactions associated with the flow (168). 



10 Conclusions 



This paper presents results of a theoretical and numerical study of many-body 
hydrodynamic interactions in suspensions of spherical particles confined be- 
tween two parallel planar walls. Our primary results include the derivation 
of transformation relations between spherical and Cartesian basis sets of so- 
lutions of Stokes equations. The transformation formulas enable construction 
of Stokes-flow fields that satisfy appropriate boundary conditions both on the 
planar walls and on the spherical particle surfaces. Using these transforma- 
tions, we have developed an efficient numerical procedure for evaluating the 
many-body resistance matrix characterizing hydrodynamic forces acting on 
suspension particles in the two-wall geometry. 



The basis sets of Stokes flows that are employed in our analysis are closely 
related to the spherical solutions introduced by Lamb and Cartesian so- 
lutions introduced by Faxen [5(j (See also section 7.4 of 5l|). By a careful 
choice of the defining properties, however, we have achieved a symmetric ma- 
trix formulation of the hydrodynamic-interactions problem. The underlying 
symmetries of the basis sets include the curl expressions linking the basis 
fields of different tensorial character, and the diagonal representations of the 
Oseen tensor in the spherical and Cartesian bases. Exploring the symmetry 
relations in our canonical formulation, the problem has been reduced to a set 
of simple explicit expressions. 



The results of our theoretical analysis were implemented numerically in an 
algorithm for evaluating the many-particle resistance matrix in the two-wall 
system. As a whole, the algorithm is quite complex, because it involves a 
large number of components. These components include constructing matrix 
elements of the Green function in terms of lateral Fourier integrals, using 
subtraction techniques to improve convergence of the integrals for configura- 
tions with widely separated particles, solving a linear system of equations for 
induced-force multipoles, and correcting the solution for slowly convergent lu- 
brication contributions. All the elements in the procedure, however, are either 
given explicitly or in terms of simple quadratures. 
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Our numerical algorithm has been used to evaluate the hydrodynamic resis- 
tance matrix for a single particle, a pair of particles, and linear arrays of 
particles confined between two planar walls. The results for the linear arrays 
indicate that the far-field flow in many-particle systems may produce signif- 
icant collective effects. A characteristic example is the large hydrodynamic 
resistance for the transverse motion of an elongated array in a narrow space 
between the walls. A simple superposition approximation in which the flow 
scattered from the walls is represented as a combination of two single-wall 
contributions fails to describe such collective phenomena. 

Our current implementation of the Stokesian-dynamics algorithm for suspen- 
sions confined between two parallel planar walls allows evaluation of hydrody- 
namic interactions in a system of about a hundred particles. For a given num- 
ber of particles, the numerical cost of the method increases with the particle 
separation (especially for p > 20H). This increase results from the oscillatory 
character of the integrands in the Fourier representation of the matrix elements 
of the Oseen integral operator. The limitation can be removed by subtracting 
several terms of the multiple-image sequence for the flow produced by the 
force multipoles induced on the particles. The subtracted contributions can 
be evaluated explicitly |27| . 

An alternative and more efficient approach is to use asymptotic expressions 
for the far-field form of the flow field produced by the force multipoles in 
the space between the walls. We have recently derived a complete set of such 
expressions. An important advantage of this approach is its simplicity — the 
asymptotic multipolar flow fields can be obtained from the solution of the 
two-dimensional Laplace's equation for the pressure field in the Hele-Shaw 
approximation. In particular, algorithms based on this method can be rela- 
tively easy generalized for periodic systems. Moreover, the efficiency of such 
algorithms can be substantially improved by applying the acceleration meth- 
ods that have been developed for Laplace's equation |52j|. We will describe 
these results in forthcoming publications. 
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was supported by NASA grant NAG3-2704 and in part by KBN grant No. 
5T07C 035 22. J. B. was supported in part by NSF grant |cts^s0348 175 and in 
part by Hellman Foundation. 



A Spherical basis 

In this appendix we list expressions for the reciprocal basis sets (16), (17) and 
(26), (27) in terms of the normalized vector spherical harmonics, as defined 
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by Edmonds 0, 



Y„_ lm (r) = af V-' +1 V 



r l Y lm (v) 



(A. la) 



Y H+lm (r) = A~V +2 V 



Y« m (r) =7,'rx V s ^ m (r). 



(A.lb) 
(A.lc) 



Here 



YU*) = n^(-l) m Pr(cos9)j m * (A.2) 
are the normalized scalar spherical harmonics, and the normalization coeffi- 
cients are 

(A.3a) 



a, = [/(2i + l)] 1/a , 

A = [(/ + i)(2/ + i)] 1 / 2 , 

7i = -i[/(/ + l)] 1/2 , 



and 



nir, 



4?r (l + m)\ 



1/2 



2l + l(l-m)\ 

The vector spherical harmonics (A.l) obey the orthogonality relations 

(5(a)Y lnm | Yi' n i m i) = Su'd nn r5 mm i . 



(A.3b) 
(A.3c) 

(A.4) 
(A.5) 



The angular functions V^ CT and W^ CT in equations (16), (17) and (26), (27) 
have the following spherical-harmonics expansions 



V£„ r = i; Y »-i+^mV ± (/;a / |a) 



(A.6a) 



W 



Ima 



(A.6b) 



The explicit expressions for the matrices V 1 * 1 at the right side of (A. 6) are 



and 



v + (0 



oli 
T+ill 



' 



2(2l+l) ai 





v-(0 



i 



21 + 1 



(Z+l)(2i+l)(2i+3)^ 



l(2l-V)(2l+l) ai ® ^ 





2(2Z+1) 



(A.7) 



(A.8) 
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Due to orthogonality relations (23) and (A. 5), the matrices W and V 1 * 1 satisfy 
the corresponding orthogonality condition 



[W ± ] t = [V*!- 1 



(A.9) 



which yields 



W+(Z) 







i 







-i 








(2+1X^+3)^-1 q (*+l)(2l+l)(2Z+3) p-1 



21 



(A.10) 



V\T(/) = (2/ + 1) 



l(2l-l)(2l+l) -1 
i+1 "i 






2(1+1) a ' 





-ihf 1 o 
o A" 1 



(A.ll) 



In the original publication [3^ and in following papers 0, H3, the basis 
functions v^ CT and w lnU7 were normalized differently. The relation between 
the spherical basis fields v^^ FS ' ) and , in the original normalization of 

Cichocki et al. and the basis introduced in the present paper is 



vw(r) = A£\>rif S) (r), v+ m<J (r) = iV^v^f s) (r), (A.12a) 

w w ( r ) = N lc ni m rw^ FS) (r) , (r) = A^% m rw^ FS) (r) , (A. 12b) 

where 

N l0 = l, N n = -(l + l)-\ /[(/ + l)(2/ + l)(2/ + 3)]- 1 . (A.13) 
B Transformation vectors X* and X 1 



The transformation vectors X*(m) and X r (m), m = —1, 0, 1 in relations (142) 
are obtained by inserting the multipolar expansion (28) into the definitions 
(140) and (141). The resulting expressions are evaluated using formulas (27), 
(A.6b), and (A. 10), which yield 



xv-r 



;i-) 1/2 



— 1 



X*(0) 



'W 2 






y/2 



X»(l) 



-1 
— i 


(B.l) 
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and 

X r (m) = -2iX*(m), m = -1,0,1. (B.2) 



C Elements of transformation matrices Tgc" an d "'"cs 



Due to the symmetric formulation of the problem, all four transformation 
matrices (90) depend on the same set of coefficients (91). Thus, it is sufficient 
to derive the explicit expression for only one of these matrices. Here we focus 
on the transformation relation (79) between the Cartesian and spherical basis 
fields and vjj^. 

To find the required transformation formula, we expand vj^ in powers of the 
radial coordinate r. The expansion can be represented in the form 

oo 

v^(r) = Ev k + i n) (r), (C.l) 

n=l 

where the fields v£^(r) are homogeneous functions of order n in r. We note 
that each term in the expansion (C.l) is itself a Stokes flow, because the linear 
operators in the Stokes equations do not couple terms with different powers 
of r. It follows that the consecutive expansion terms can be represented as 
combinations of the spherical basis solutions (17) with / + a — 1 = n. For the 
pressure solution v^ 2 , this representation can be expressed as 

v +(«) _ „(«) i „(") i ,,( n ) (n 9\ 

V k2 — U k2 + U kl + U k0 > y^-^) 

where 

uL +CT_1) = E 4>L, (c.3) 

m=—l 

Comparing the above expressions with equations (79), (89), and (90a) yields 



a 



kcr 



i m (27ck)- 1/2 k l+a+1 e- im ^a CT , (C.4) 



where the coefficients a a correspond to the coefficients a, b, c in equation (90a), 

a = c, a± = 2b, a<i = 4a. (C.5) 

The remaining components v^™- 1 and v^"- 1 can be related to the flow fields 
(C.3) by applying the curl operator to both sides of equation (C.l) with a = 2 
and a = 1. After inserting decomposition (C.2) and using curl relations (57) 
and (206), we collect terms corresponding to the same power of r, which yields 

v^-^li^Vx^g + uff), (C.6a) 
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v^ 2) = -^- 2 Vx(Vx U S). (C.6b) 

The above results are consistent with the triangular structure of the trans- 
formation matrix (90). To evaluate the coefficients a' fc ™, we reduce expressions 
(C.2) and (C.6) to equivalent relations between appropriately defined har- 
monic scalar fields. For the spherical basis flows Vi ma we have 

viL = ilA, (c.7) 

where 

<P+ m (r)=r l Y lm (i), (C.8) 
and the operators Lf a are given by 

Lfo = V, (C.9a) 

Lfi = i(/ + l)- 1 rxV, (C.9b) 

Lf 2 = [(i + 1)(2Z + ^[-h- + \{l + 3)r 2 V]. (C.9c) 
The analogous expressions for the Cartesian basis fields are 

vt = il$t, (cio) 

where 

<2>+(r) = (32 7 r 2 A;)~ 1 / 2 e ik -^. (C.ll) 
The operators L£ and Lj^ are given by the expressions 

Lfeo = k-'V, (C.12a) 



Lg = 2ik~ 1 e z x'V, (C.12b) 
and the operator Lg is given by 

I& = Lg + Lg, (C.12c) 

where 

Lg = fc-iv, Lg = -2e, + 2zV. (C.12d) 
The above relations can easily be verified using expressions (17) and (A. 7) for 
the spherical basis fields and expressions (51) for the Cartesian basis. 

The radial-expansion components in the decomposition (C.l) of the 

Cartesian basis flows can be determined by applying the operators Lj^ to the 
expansion of the Cartesian scalar field (C.ll) in powers of r 

oo 

*i = E*£ in) , ( c - 13 ) 

n=0 
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where 

< {n) (r) = (32 7 r 2 A;)- 1 / 2 (ik ' P J" MTt . (C.14) 

Inserting relations (C.10) and (C.13) into the expansion (C.l) and collecting 
terms corresponding to a given power of r we find 

v k + (n) = lW 1 ' v k + W = t C kl K (n+1 \ (C15a, b) 

and 

v k + 2 W = t%< in+1) + Lg<<">. (C15c) 
With the help of relations (C.7)-(C.9) for the spherical basis fields, the flow 
fields in equations (C.6) can be represented in a similar manner, 

uL +CT_1) = LfX. ( °, ^ = 0,1,2, (C.16) 

where 

= E 4Wn, (C17) 
m=— 2 

according to equation (C.3). 

A closed set of equations for the scalar functions is obtained by inserting 
relation (C.16) into (C.2) and (C.6), using (C.15), and employing the curl 
identities 

iVxLf 2 !Pi = Lf^, iVxLf^ = tf W t (C.18a, b), 

where is an arbitrary solid harmonic of the order /. The above expressions 
correspond to the curl identities (19) for the spherical basis fields, and can 
be verified using relations (C.9). The equations for the scalar functions 
derived by this procedure are 

Lf ^=4A; 2 LK«, (C.19a) 
L^+« = 2*L°*J (,) - tlnKt 1 ', (C19b) 

W = £g*f (,) + - Lf.n^'- 1 ) - Lf_ 22 </- 2 ). (0.19c) 

The above equations can be explicitly solved for the unknown fields ■ Us- 
ing expressions (C.9) and (C.12) for the operators Lp CT and L^ CT , and simplifying 
the results using relation (C.14) for the field we find 

#+ (,) =4fc*£ (,) , (C.20a) 



47 




(C.20b) 



and 



+(0 fe[(2Z 2 - 4/ + 3)(Lr + + 2(1 - 2)z(ia; + g) - 2(1 -!)(/- 2)y 2 ] +a _ 2) 
k0 " Z(Z — 1)(2Z — 1) k 

(C.20c) 



where k = ke x is assumed. 

In the final step of our derivation, we recall that the functions are solid. 

harmonics of order /, according to expressions (C.8) and (C.17). To obtain the 
expansion coefficients a 1 ^ in relation (C.17) for even values of the parameter 
I + m, we evaluate both sides of (C.20) on the plane z = and compare 
the coefficients of the angular Fourier modes e im< t>. In the case of odd values 
of the parameter / + m, a similar analysis is performed for the derivative of 
both sides of equations (C.20) with respect to the coordinate z. The analysis 
yields the quantities a 1 ^ in the form (C.4), with the coefficients (C.5) given 
by expressions (91). 



D Large k behavior of integrands 5^/(k) 

In this appendix we derive the asymptotic expression (153) for the large k 
behavior of integrand 5^/(k). According to equations (133) and (149), the 
decomposition (150) of the integrand ty(k) corresponds to the separation 



of the two-wall scattering matrix (136) into the 0(1) diagonal contribution 



Z TW (k) = Z + 5Z TW (k) 



(D.l) 



(D.2) 



Z 



■w 



and the correction of the form 



5Z TW (k) 



Zq • 5tw(^) • Z-£w(k) 



(D.3) 



where 







S++(-kH) 



S TW (k) 



(DA) 



S c -(kH) 
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Taking into account relation (70) we find that 

5Z TW (k) ~ e~ kH , (D.5) 

which implies that 

5Z TW (k) < Z , k > 1. (D.6) 

Inserting the decomposition (D.l) into equations (133) and (149) yields 

6*(k) = 5g TW (k; Ima \ I'm' a')k l+v+(7+(T ' ~ 2 J m ,^ m (k Pl] ), (D.7) 

where 

SgTw(k;lm \ I'm') = -t sc (lm) -S iW (k) -6Z TW (k) -S Wj (k) ■ t C s(l'm'). (D.8) 

The asymptotic expression (153) is obtained by inserting relation (D.3) with 
Z TW (/c) ~ Z into (D.8), and using equations (68), (70) and (127). Evaluation 
of the slowest-decaying term yields (153) with 

Aij = mm(Z iL + Zuj, Z jL + Z m ) + H, (D.9) 

which is equivalent to (154). 

We note that the convergence of the integrand (149) can further be improved 
by subtracting from the two-wall scattering matrix Zpw several terms in the 
expansion 

oo 

ZTw(k) = £(-irZ • [StwW • ~ZoY- (D.10) 

One can show that the subtracted terms correspond to consecutive reflections 
of the flow field from the walls. Thus, these terms can be evaluated without 
numerical integration using the image-representation formulas derived by two 
of us Q. 
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